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I. INTRODUCTION 

In the last years there has been significant progress in the 
analysis and characterization of the dynamics of processes un- 
derlying the time series of complex dynamical systems 0}[3]. 

If the temporal evolution of a quantity X t can be described 
by a Langevin equation, it is possible to extract drift and dif- 
fusion functions of the underlying stochastic process from a 
given time series. This can be done because the moments of 
the conditional probability densities of X t + T \x t =x can be re- 
lated to these functions. 

Since this approach was introduced ||3}{8l it has been suc- 
cessfully carried out in a broad range of fields. For example 
for data from financial markets [9], traffic flow ifTOl . chaotic 
electrical circuits lfTTl[T2l . human heartbeat [13], climate in- 
dices lfl4l [TBI , turbulent fluid dynamics fl6l . and for elec- 
troencephalographic data from epilepsy patients ifFTl [181 . 

Real-world data, however, also give rise to some problems. 
One of them is, that experimental data are only given with a 
finite sampling rate. So methods had to be proposed to deal 
with the effects arising from this fact |[T9Tl2Tl . 

Another problem is the virtually unavoidable measurement 
noise [ 1., 2 1-23 1. In the presence of measurement noise Y t the 
values of X t or any of its probability densities are no longer 
accessible, but only X£ = X t + Y t and its density distribu- 
tions. 

Recently an approach has been presented which allows the 
estimation of drift and diffusion functions in the presence of 
strong delta-correlated, Gaussian noise [24, 25 1. Starting with 
initial estimates for the noise strength and the drift and dif- 
fusion functions a functional of these unknowns is iteratively 
minimized. 

The aim of this paper is to introduce an alternative approach 
to the treatment of strong measurement noise. This approach 
is able to deal also with exponentially correlated, Gaussian 
noise. The basic idea is, not to look at the conditional mo- 
ments in the first place but at the joint probability density 
p(xi,x 2 ,t) of pairs (X t , X t+T ). If the measurement noise 
is independent of X t , then (X t , X t + T ) and (Y t , Y t + T ) are in- 
dependent random variables and the joint probability density 
p*(xi,X2, t) of their sum , X% +T ) is given by the convo- 
lution of p and py, where py (xi,X2, r) is the joint probability 
density of (Y t ,Y t+T ). 

The noise is assumed to be Gaussian and the Gauss func- 
tion has some special properties with regard to convolution 



and Fourier transform. This makes it possible to extract the 
noise parameters from the moments of p*. Furthermore the 
abovementioned relation between the conditional moments 
and the unknown functions can be transformed into a relation 
in Fourier space. This allows polynomial approximations of 
the drift and diffusion functions to be extracted using purely 
algebraic relations between quantities that can be calculated 
directly from a given, noisy time series. 

This paper is organized as follows: Section [II] is devoted 
to the noise-free stochastic process, the definition of its joint 
probability density and expressions for the moments of this 
density in terms of a Taylor-Ito expansion. Section III pro- 



vides the properties of the measurement noise under consider- 



ation and in section IV expressions for the moments of a noisy 
process will be derived. In section|V]these expressions will be 
used to extract the parameters of the measurement noise and 
in section [VT] to extract polynomial approximations for drift 



and diffusion functions. Finally in section VII the results are 



applied to some synthetic time series. The used properties of 
the Gauss function and further computational details are given 
in appendices [A]to |C| 



II. STOCHASTIC PROCESS 

Let X t be a stochastic process that can be described by a 
time-independent Ito Langevin equation 

dX t = D^{X t )dt + yjDW(X t )dW t , (1) 

where D' l > and are the Kramers-Moyal coefficients of 
the corresponding Fokker-Planck equation and dW is the in- 
crement of a standard Wiener process (< dWdW >= dt). 

The above equation cannot be solved analytically in gen- 
eral. To get finite-time step approximations for the temporal 
evolution of X t , a Taylor-Ito expansion can be applied. 



X t+T =X t + h(X t ,r) 



R 



(2) 



The function h contains the expansion terms and R is a re- 
mainder that will be dropped when using the approxima- 
tion. The above equation allows the numerical integration of 
Eq. ([T]i but will not be needed in the following. What will be 
needed are approximations for the first and second moment of 
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X. 



t + T 



x. 



mi(xi,r) = h 1 (x 1 ,T)m (x 1 ) 



(lib) 



< (X t+T - X t )\x t =x > = hi(x,r) + Ri (3a) 
< (X t+T - X t ) 2 \ Xt=x > = h 2 (x,r) + R 2 (3b) 

A Taylor-Ito expansion also provides expressions for these 
moments. For an expansion of order k the functions hi and 
h 2 are polynomials of order k in r (without a constant part) 
and the remainders Ri and R 2 are of order r fe+1 . 



hi(x,r) = ci(x)t + . . . + c k (x)r k 
h 2 (x,r) = d\{x)T + . . . + d k {x)r k 



(4a) 
(4b) 



The polynomial coefficients generally are functions of 
and D 1 - 2 ^ and their derivatives. The first order approximations 
that will be used later are given by 



hi(x,r) = D {1) {x)t 
h 2 {x,r) = D (2) (x)t. 



(5a) 
(5b) 



A detailed description of the Taylor-Ito expansion and its mo- 
ments can be found in ||26l . 

The probability density of X t and the joint probability den- 
sity of (X t ,X t + T ) are formally given by 



lim — 

T-)-oo T 



p(x 1 ,x 2 ,T) = lim - 

1 — >oo 1 



S(X t - xt)dt 



5{X t - Xl ) 
8(X t+T - x 2 )dt. 



Let the conditioned moments nij be defined as 

/oo 
(x 2 - x 1 yp(x 1 ,x 2 ,T)dx 2 . 
-oo 



(6a) 



(6b) 



(7) 



(These moments are not to be confused with the conditional 
moments which usually are defined as nij/m ). Inserting 



Eq. ( 6b I the moments can be expressed as 

1 f T 

nij{xi,T) = lim — / 5{X t -X\) 

x (X t+T - x x ) 3 dt. 



(8) 



Using 



lim i [ 6(X t - Xl )f(X t )dt = 

< f(X t )\ Xt=Xl > p(xx) 



(9) 



leads to 

mj(ii,r) 



< {X, 



t+T 



x t y\ Xt=xl > p(xi) (io) 



Inserting Eqs. Q and omitting the remainders the final ap- 
proximations of the moments up to order two are therefore 
given by 



m (xi) = p(xi) 



(11a) 



m 2 (xi,r) = h 2 {xi,r) m (xi) 



(11c) 



These are the equations that establish the relation between the 
observable quantities nij and the unknown coefficients con- 
tained in hi and h 2 within the precision of the chosen Taylor- 
Ito expansion. 



In the absence of noise Eqs. (Hi can directly be used to 
determine and D^ 2 \ If e.g. the above mentioned first 
order scheme was chosen, the equations read 

^4 = r^)(x l) + 0(r 2 ) 
m (xi) 

^4 = tD( 2 Hx,) + 0(t 2 ). 
m (xi) 



III. MEASUREMENT NOISE 

Let Yt be exponentially correlated Gaussian noise as pro- 
duced by an Ornstein-Uhlenbeck process. 



dY t = -aY t dt + bdWt 



(12) 



This equation can be solved analytically (see e.g. Il27l ). The 
probability densities of Y t can be expressed either in terms of 
a and b or, equivalently, in terms of the 'macroscopic' param- 
eters T (characteristic time scale) and a 2 (variance) with 



1 a b 2 

T = - and a 2 = —. 
a 2a 



(13) 



The unconditioned distribution of Y t will be denoted by K 
because it will serve as a convolution kernel later. It is given 
by 



K(x) 



1 



V2 



(14) 



7T(J 



The joint probability density of (Y t ,Y t+T ) can then be ex- 
pressed as 

Py(x!,x 2 ,t) = K( Xl )-j=e n • ; ,(15) 

where the decay-function 

M(r) = e"* (16) 
and the auxiliary quantity 

s 2 = a 2 (l-/i 2 (r)) (17) 



have been introduced for notational simplicity. The parameter 
T will also be called 'the' correlation time. This is motivated 
by 



< Y t Y t+T > 



(18) 
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IV. NOISY STOCHASTIC PROCESS 



+2(1 - p)a 2 K * m 



(22c) 



Let X* — X t + Y t be the sum of a stochastic signal X t 
and measurement noise Y t as introduced in sections [TT] and 
[Hi] respectively. Then, because X and Y are independent, 
the joint probability density p* of the pairs of noisy variables, 
(X% , X t * +T ), is given by the convolution of p and py . 

p*(x 1 ,x 2 ,r) = py(x 1 ,x 2 ,t) * p(x 1 ,x 2 ,r) 

Py{x\ - x[,x 2 - x' 2 ,t) 

x p(x' 1 ,x' 2 ,T)dx' 1 dx' 2 (19) 

Instead of the conditioned moments rrtj only their noisy coun- 
terparts to* can be determined. 



oo />oo 



OO — OO 



TO* (xi,r) 



(x 2 - xi) 3 p*(xi,x 2l T)dx 2 (20) 



These can be expressed as follows (see appendix [B}. The in- 
dex of the variable xi and the function argument of p(r) are 
omitted for notational simplicity. 

TOq(ie) = K(x)*mo(x) (21a) 

toJ(x,t) = K(x) * mi(x, t) 

-(1 -fi)(xK(x)) *m (x) (21b) 

m 2 (x,T) = K{x)*m 2 (x,T) 

—2(1 — fi)(xK(xj) * mi(x, t) 
+ (1 -p) 2 (x 2 K(x)) *m (x) 

+s 2 K(x) * mo(x). (21c) 

These relations have already be derived earlier for the case 
of delta-correlated noise and have been formulated with con- 
ditional probabilities [24 J. The above convolutional notation 
makes it more obvious how to proceed further. Because K is 
a Gauss function with variance a 2 , the terms x 3 K{x) can be 
expressed as (see appendix|A|): 

xK(x) = -a 2 d x K(x) 
x 2 K{x) = (a 2 + a 4 dl)K(x) 



Because of (d x f) * g — d x (f * g) this leads to (using Eq. ( 17 1 
and omitting arguments now completely): 



to,q = K * too 

toJ = K * mi + (1 — n)p 2 d x (K * mo) 

m 2 — K * to 2 + 2(1 — fi)a 2 d x (K * toi) 
+ (1- n) 2 a A d 2 x {K*m ) 



(22a) 



(22b) 



Substituting K * too and K * toi and using Eqs. (22ai and 
( |22b| yields 



TOq = K * too 

toJ = K * toi + (1 — p)a 2 d x m*Q 

m* 2 = K * m 2 + 2(1 — u,)a 2 d x m* l 

-(1- [i) 2 a 4 d 2 m* +2(1- p)a 2 m*. 



Using Eqs. ( 1 1 1 to express m\ and to 2 finally gives 

TOq = K * too 

to^ = A" * (/iimo) + (1 — p)a 2 d x rriQ 

m 2 = K * (h 2 mo) + 2(1 — p)a 2 d x m\ 

-{l-tfa i %m* Q + 2(\- t i)<T 2 n%. 



(23a) 
(23b) 

(23c) 

(24a) 
(24b) 

(24c) 



Eq. (24b i can be used to extract the parameters of the mea- 



surement noise without determining the Kramers-Moyal coef- 
ficients £>W and D^ 2 \ This will be done in the next section. 



V. EXTRACTING MEASUREMENT NOISE PARAMETERS 



Multiplying Eq. (24b i by a weight function ty(x) and sub- 



sequently performing an integration with respect to x leads 
to 



A(r) 



A(r) 



B = 



(K(x) * (hi(x, T)m (x)) 1 $(x)dx 

>o 

+(1 - p(r))a 2 B (25) 
m* (x, r)^f(x)dx (26) 

X) 

ml(x)(d x ^(x))dx (27) 



In the last step an integration by parts has been applied, as- 
suming \l/77ig — > when |x| — > oo. Using Eq. (4ai the re- 



with 



and 



maining integral in Eq. ( 25 1 can be expressed as a polynomial 



{K(x) * (hx(x, T)m (x)}^(x)dx 

fe /*°° 

/ {if (as) * (cj(a:)mo(a:))}*(a:)da; 
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3=1 



Finally, after a division by B, Eq. (25 i reads 

z(t) = (1 - p(r))a 2 + Cir + C 2 t 2 + 

with 



z(r) 



A(r) 
B 



and 



(28) 



(29) 



(30) 



Once a weight function has been chosen, z(t) can be de- 
rived from experimental data. Eq. (|29]> then can be used to fit 



the unknown parameters a , T (contained in p = e~T) and 
Cj on the right-hand side. This can be done for example by 
an iterated least-square fit. 

The simplest choice if? = x is not the most accurate one, 
especially in the presence of heavy tails. It puts a high weight 
on the tails of m\, where density is low and fluctuations of 
experimental data thus high. Choosing 9 X \& to be a rough, 
piecewise steady approximation of the density TOq overcomes 
this problem and also assures B ^ 0. 

No binning has to be applied in the evaluation of z(r), if 
the values of X* are given at constant time-increments At. 
Let Xi = X*(ti) with U = to + iAt. Then, instead of us- 
ing box-counting techniques to approximate the density p* , 
the temporal evolution of X% itself can be approximated by a 
piecewise constant function: 



x; = 



(ti<t< ti+i) 



(31) 



The density distribution of , X^ +Tk ) with r k = k At can 
be evaluated using Eq. ( |6b| and is given by a sum of Dirac- 
distributions (the coordinates x\ and x' 2 are used to avoid con- 
fusion with the values xi of the given time series). 

1 N 

p*(x' 1 ,x' 2 ,T k ) = —^2S(x l -x' 1 )S{x l+k ~x / 2 )(32) 
»=i 

The moments of p* evaluate to 

1 N 

m^x'^Tk) = —^2{x i+k -x' 1 ) : >8{xi-x' 1 ) (33) 



what results in 



Mjk) = y^S X i+k - Xj)^(Xj) 

i=l 
1 N 

13 = -vE*^) 



where denotes the derivative with respect to x. The values 
z(rfc) derived from a given time series Xi can therefore be 
expressed as 



z(r k ) 



(35) 



To be able to successfully perform a fit it must be possible to 
distinguish between the polynomial CjT^ and the function 



1 — within the given range of increments r. This means 
that a good polynomial approximation for e~T should require 
a higher order than the polynomial defined by the Cj . The pro- 
posed method is therefore limited to measurement noise with 
a correlation time T considerably smaller than the time scale 
of the underlying stochastic process. An illustrative example 
is given in Fig. [T] 
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FIG. 1: Values z derived from experimental data (Ornstein- 
Uhlenbeck process) with and without measurement noise (a = 1 
and (7 = respectively). The correlation time T is given in units of 
At. The process time scale is lOOAt. 



VI. EXTRACTING DRIFT- AND DIFFUSION 
COEFFICIENTS 



Going back to Eqs. ( |24| i and putting the convolution terms 

(36a) 



on the right-hand side gives 

ml — (1 — fi)a d x rriQ = K * (himo) 

m* 2 — 2(1 — p)a 2 (d x ml + mp) 

+ (1 - n) 2 a 4 dlm* = K * (h 2 m ). (36b) 
To get rid of the derivatives of m* a Fourier transform 



F : /(*) f(u) 



c f{x)dx 



(37) 



is applied. In Fourier space Eqs. (36 1 now read 



m* — Miloijiq = 71 (38a) 

m* 2 - 2M(ml + iujm\) - M 2 cj 2 m* = 72 , (38b) 



(34) 

with the right-hand sides 



71 = F{K *(/iim )} 

72 = F{K * (h 2 m )} 

and the abbreviation 

M = (1 - /i)cr 2 . 



(39a) 
(39b) 



(40) 



If £>W and can be expressed as polynomials in x, 



Eqs. (38 1 can be used to extract the corresponding polynomial 
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coefficients. Let and be given by 

Ni 



D (2) 



E^ 

3=0 

N 2 



x J 



(41a) 



(41b) 



Choosing the first order approximation Eqs. (|5]l for hi and h 2 , 
the right-hand sides of Eqs. (|38j) read 



7i = rFjiC* (L> (1) to )} 



rV^ajF (if * (a;- 7 to )) 



(42a) 



72 = T 



F^K * (Z? (2) to )| 



iV 2 



3=0 



(42b) 



So for polynomial drift and diffusion functions the problem 
of expressing ji and 72 reduces to finding expressions for the 
transforms of K * (x^rno). It should be noted here that this 
would also be the case for higher order approximations hi and 
h 2 . 

Again because K is a Gauss function with variance a 2 , the 
Fourier transform of K * (a^m-o) can be expressed as lin- 
ear combination of the transforms of x k (K * too) (see ap- 
pendix [A}. Since K * too = TOg these transforms can be de- 
rived from the noisy data. Using the shortcuts 



<f> 3 = F(x j (K*m )) = F(x j m*) 
Fj = F (K * (Vto )) 

the relation between $ and F is given by 



F; 



k=0 



with 



i 3 J2\a jk \a 



3 + k^k 



fe=0 



and 



; fc = ~ a (j-i)o-i) + s + i) a (j-i)(/c+i) 



aoo = 1 

a jk = 



(j<0,fe<0,i<fc). 



(43) 
(44) 



(45) 



(46) 



(47) 



Eqs. (38 1 now read 



JVi 



ij — MiunfiQ = t CLjFj 
3=0 



(48a) 



to*. - 2M(m* + mjtoJ) - M 



,2 ~ * 



N 2 

•E 



(48b) 



If the noise parameters have been extracted according to sec- 
tion [v) then M = (1 — m) 

a and Fj are known and the co- 
efficients cij and fej can be extracted by a least square fit. It is 



also possible to use Eqs. (48 1 to fit noise parameters and poly- 



nomial coefficients simultaneously by an iterated least square 
fit. But this has turned out to be less accurate. 

Adding an increment dependency to the polynomial coef- 
ficients however allows some of the approximation errors to 
be absorbed by the additional parameters and has shown to 
increase accuracy. 



rb) 



(49a) 
(49b) 



With this extension the final equations for parameter fitting 
are given by 



1* — Miu>rh,Q 



-E< 

3=0 



aj +Ta))F 3 



(50a) 



m* 2 - 2M (to* + iioml) - M 2 uj 2 m* 

N 2 

= rY,%+rb))F . 

3=0 



(50b) 



To be able to actually perform a fit, a finite number of values 
for u> and r has to be chosen. It is sufficient to restrict the 
choice of lu to positive values because the above equations are 
symmetric (conjugate complex) in to. The following heuristic 
approach has shown to work quite well: 



Define an upper bound lu* by 



\m*J 2 duj = 0.99 



'0 



■duj. 



(51) 



Chose values ljj equally distributed in [0,w*]. The 
actual number of values is of minor importance but it 
should be high enough to properly sample TOq. 



VII. APPLICATION TO NUMERICAL DATA 

In order to test the accuracy of the proposed method, syn- 
thetic stochastic signals have been generated by numerical in- 
tegration. As a first test-case an Ornstein-Uhlenbeck process 
has been chosen. Its drift and diffusion functions have been 
defined as 



D^(x) = -x 
D^(x) = 2. 



(52a) 
(52b) 



The generated time series will be referred to as data set A. It 
consists of 10 6 data points at time increments At = 10~ 2 . 
The integration has been performed with a timestep size of 
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FIG. 2: Probability density function for data set A. a is the standard 
deviation of the superimposed noise. Drift and diffusion functions 
are given by D' 1 ^ = —x and I)' 2 ' = 2. 



FIG. 4: Estimated polynomial coefficients of D^ 1 ' for data set A with 
ao and a\ according to Eq. l |41[ >. 
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FIG. 3: Absolute error Act = a — a in the estimation of the mea- 
surement noise of data set A. The dashed lines indicate the expected 
lower bound of the standard deviation of Act. 

10~ 4 using the Euler scheme. The process-timescale of 1.0 
corresponds to r = 100 in units of increments of data points. 

This time series has been superimposed by additional mea- 
surement noise with standard deviation in the range a = 
0.0 . . . 2.0 and correlation time T = 0.0 (delta-correlated). 
Some of the density histograms of the resulting signals are 
shown in Fig. [2] 

The noise parameters have been extracted according to sec- 
tion[V] The maximum increment has been r = 60. The abso- 
lute error Ac in the estimation of the strength of the mea- 
surement noise is presented in Fig. [3] as a function of the 
true measurement noise a. The results are accurate within 
the achieveable precision of the finite time series. This means 
that the measurement noise alone, without any stochastic pro- 



FIG. 5: Estimated polynomial coefficient bo of D' 2 ' for data set A 
with bo according to Eq. (41) . 



cess, cannot be estimated much more accurate: For a sample 
of N data points of Gaussian, delta-correlated noise the esti- 
mation S of the true standard deviation, a, obeys a variance of 
approximately a 2 ^ (see appendix|cji. A lower bound for the 
standard deviation of Act is therefore given by a—k=. These 
limits are indicated in Fig. [3] as dashed lines. 

Next the drift and diffusion functions have been fitted 
choosing the polynomial ansatz = ao + a\x and — 
b . As maximum increment a value of r = 25 (again in units 
of increments of data points) has been chosen. The results are 
shown in Figs. |4] and [5] and are in good aggreement with the 
true values (ao = 0, a\ = —1, bo = 2). 

After this basic test-case a more general case has been cho- 
sen. It includes multiplicative noise, heavy tails and a density 
distribution which is not symmetric. Also the superposition of 



7 




FIG. 6: Probability density function of data set B. a is the standard 
deviation of the superimposed noise. Drift and diffusion functions 
are given by = 1 - x and D (2) = 2 ~ 2x + 2x 2 . 



noise with finite correlation time has been investigated in this 
setup. Drift and diffusion functions have been defined as 



DW(x) = 1-x 
D {2 \x) = 2-2x + 2 



x 2 . 



(53a) 
(53b) 
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FIG. 7: Results for data set B. Estimations of the polynomial co- 
efficients a,j and bj and absolute error Act in the estimation of the 
measurement noise. The superimposed noise is delta-correlated. 



The generated time series will be referred to as data set B. 
Again it consists of 10 6 data points at time increments At — 
10~ 2 . The integration has been performed with a timestep size 
of 2.0-10 -5 using the Euler scheme. The standard deviation of 
the data points is approximately 1.83 and the timescale of the 
process has been estimated via the autocorrelation function to 
be about r = 100 in units of increments of data points. 

In a first step again delta-correlated measurement noise 
with standard deviation in the range a = 0.0 ... 2.0 has been 
superimposed to the process. Some of the density histograms 
of the resulting signals are shown in Fig. [6] 

As maximum increment for noise-fitting r = 60 and for 
coefficient-fitting t — 25 has been chosen. The polynomial 
ansatz for drift and diffusion functions is given by £>W = 
a + <X\X and = b + b\X + b 2 x 2 . The results are sum- 
marized in Fig. [7] 

In a second step measurement noise with a correlation time 
T = 2.0 (in units of data point increments) has been super- 
imposed. The results are given in Fig. [8] and obey higher 
fluctuations as in the delta-correlated case. The correlation 
time of the measurement noise is estimated nicely except for 
small amplitudes of a. To check the assumption that this is 
caused by the finite-sample fluctuations of z(t) a larger data 
set, Bl, has been generated. It consists of 10 7 data points and 
its analysis should show less variability in all estimated quan- 
tities. This is indeed the case. The results are shown in Figs. 
|9]and[l0] 



VIII. CONCLUSIONS 

A new procedure has been described to analyze stochas- 
tic time series that are superimposed by strong measurement 
noise. The algorithm is able to cope with exponentially cor- 
related noise and accurately extracts strength and correlation 
time of the measurement noise as well as the parameters defin- 
ing the drift and diffusion functions of the underlying stochas- 
tic process. This has been shown by the analysis of synthet- 
ically generated time series. The chosen stochastic processes 
include the cases of multiplicative noise and heavy tailed den- 
sity distribution. 

The computational costs of the algorithm are quite low. It 
takes less than a minute to analyze a data set of size 10 7 data 
points on a usual desktop PC. 

A first order Taylor-Ito expansion for the moments of the 
finite time differences of the process variable is used in the 
current implementation. It is straightforward to extend the 
algorithm to also take higher order terms into account. This 
should extend the range of time increments that can be used 
for the analysis and thus should increase the accuracy of the 
results. This has to be done in the future. 

Another future task is the extension to higher dimensional 
processes. Up to now it is an open question if a general mul- 
tidimensional Omstein-Uhlenbeck process as a model for the 
measurement noise can be used or if some restrictions need to 
be imposed on it. 
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FIG. 9: Results for data set Bl. Estimations of the polynomial co- 
efficients a,j and bj and absolute error Act in the estimation of the 
measurement noise. The superimposed noise is delta-correlated. 



choice. 



FIG. 8: Results for data set B. Estimations of the polynomial coef- 
ficients cij and bj, correlation time T and absolute error Act in the 
estimation of the measurement noise. The superimposed noise has a 
correlation time of T = 2.0. 
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Appendix A: Gauss functions 

1. Definitions and basic properties 

Let Fourier transform and convolution be defined as below. 
For notational simplicity the 'hat' syntax will be used to de- 
note the transform of single functions. For more complex ex- 
pressions the functional form F(. . .) will usually be the better 



/-t-oo 
e-^ x f(x)dx (Al) 
-OO 

/+oo 
f(x - x')g(x')dx (A2) 
-OO 



Above definitions imply the following properties: 



F(dif(x)) 




F(xif(x)) 


= W/M 


F(fC-)) 


= \a\f{au) 


F(f(x)*g(x)) 




dl(f(x)*g(x)) 


= (dif(x))*g{x) 




= f{x) * (dig(x)) 



(A3a) 
(A3b) 
(A3c) 
(A3d) 

(A3e) 
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The values of the polynomial coefficients can be determined 
by a recursion formula: 

Ojfc = — ay_i)(fc_i) + (fc + l)oy_i)(fe + i) 

a 00 = 1 (A6) 

Ojfe = (j < 0, fc < 0, j < k) 

It can also be shown, that the coefficients are zero, if (J + 
k) or (j — k ) is odd. Otherwise, within the restrictions of 
formula ( A6 1, djk is non-zero and its sign is given by P +k . 



It will be useful later, to express the derivatives of ^ in a 
similar way. The proceeding is the same as for G and leads to 
the recursion formula dA8b. 



1 



G(x) 



.fc=0 



1 



G(x) 



(A7) 



Ojfc = 5(j-_i)(fc_i) + (fc + l)o( 3 -_ 1 )( fc+1 ) 
aoo = 1 

ajfe = (j < 0, fc < 0, j < k) 



(A8) 



FIG. 10: Results for data set Bl. Estimations of the polynomial 
coefficients a,j and bj, correlation time T and absolute error Aa in 
the estimation of the measurement noise. The superimposed noise 
has a correlation time of T = 2.0. 



2. Derivatives and polynomial products of Gauss functions 



Let the Gauss function G be defined as below. 



G(x) = ^Le^* 2 



G(u) = e 



= e 2- = 



^2ttG(u}) 



(A4a) 
(A4b) 



All derivatives d 3 x G(x) have the form Pj(x)G(x), where Pj 
is a polynomial of order j in x. 

d x G(x) = {-x}G{x) 

d 2 x G(x) = {-l + x 2 }G(x) 



diG(x) = \J2a jk x k \G(x) 



(A5) 



, k=0 



Having a closer look at the coefficients, it turns out, that a,jk 
and cijk only differ in sign. 



a.jk — \o>ik\ 



(A9) 



To get explicit expressions for x^G(x) in terms of derivatives 
of G, a Fourier transformation is applied to Eq. ( A5 I. 



{iu)) j G{uj) = \J2a jk (id u ) k \G(u 



(A10) 



>fc=0 



Using Eq. ( A4b i and substituting oj by x, finally yields: 

x j G(x) = lj2i k - j a jk d*\G(x). (All) 



.fc=0 



3. Scaled Gauss functions 



This can also be expressed for scaled Gauss functions hav- 
ing standard deviation u. The Fourier transform of such a 
function is given by 



F(G(-)) = aG(au). 
a 



(A12) 
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The relations between derivatives and polynomial products 
then read: 



dlG(-) = ^a^'-V G(-) (A13) 
! ' U=o J ' 



4 m = {t ]ajkla ~ j ~ kxk }^ 



(A14) 



iG & = E^-W + ^ \G(l) (A15) 



. fe=0 



4. Convolutions 

The above properties can be used to express convolution s of th e for m (x 3 G) * f and G * (x 3 /) in ter ms of t he 'raw' convolution 
G * f and its polynomial products xP{G * f). Using Eqs. (A3ei and (A15 i immediately leads to Eq. (A16l. 



(x'GU)*/(*) = NE iH ^^| G ( ff )J */(*) 



To get an expression for G * (x 3 f) it is useful to look at the term G ■ d 3 f first. 



G(-)d 3 J(x) = G(-)dl\-±-G(-)f(x) 

1 



■'■ ^ / 3 

k=0 ^ 



a; X-*/ x \ 1 x 



ni~ k ^ \d*{G{^)f{x)} 

3 / \ ( i—k 



G(f) 



G OE J Ei«H*-*)«k- 

fc=0 ^ ' u=o 
J / -\ j-k 



E L EW" 



-i+ k - l x l 



k=0 v 7 Z=0 



^jGT|)^{ G ^^} 
^{G(^)/(x)} 



This allows the transform of G * {x 3 f) to be written as follows: 

F (G(^) * {x 3 f{x))) = <7£(<7w)(t0 u )*/(w) - iVV^ • G(tru)d>J(u) 

= iV^E (fc) E l^-^k^^ {g(™)/»} 

fc=0 ^ ' 1=0 

= E (3E rl i%-^i ffr ^' F (^( G (^/w)) 



/c=0 x ' (=0 

Omitting arguments and introducing the auxiliary quantity <pj this can be expressed as 



k=0 



F (G * (aJ/)) - E 1 ^-* F ( x "i G * /)) . ^ = f E 



(A16) 



(A17) 



(A18) 



fe=0 
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Appendix B: Conditioned moments m* 



After inserting Eqs. ( 15 i and ( 19 1 and interchanging the order of integration Eq. (20 1 reads 

/CO />oo 
K(x\ — x' 1 ) / p(x' 1 ,X2,T)Ijdx' 2 dx' 1 , 
-co J — CO 

where 



Ij = 



Substituting 



(x 2 - Xi) 3 



V2ns 2 



1 f X2-f-( T ){ x l- x \)- x 2 \ 



dxo 



%2 - M( T )( a; l _ x 'l) - x 2 

Z = 

s 

gives Ij — f^oisz + (x' 2 — x[) — (1 — h(t))(xi — x^jy ^ dz and therefore 

1 



h 

h 



(x' 2 - x[) - (1 - /i(T))(a;i - xi) 

s 2 + (4 - xi) 2 - 2(1 - M (t))(xi - 4)(^2 " A) + (1 ~ M(^)) 2 (^i " ^i) 2 
Now Ij = p(x' 1 ,x' 2 ,T)Ijdx' 2 can be evaluated using the definition of mj (Eq. (|8j). 

Iq = m (xi) 

7( = mi(xi,r) - (1 - /i(r))(a;i - x[)m (x[) 

I' 2 = m 2 (x' 1 ,r) - 2(1 - h(t))( Xi - z'JmiK, r) + {(1 - mW) 2 (^i - ^i) 2 + s 2 } m (x\) 
Using (x — x') k f(x — x')g(x')dx' = (x k f(x)) * g(x) the noisy conditioned moments read 

m* a (xi) = K(xi) *m {xi) 

ml{xi,T) = if(xi) *mi(a;i,r) - (1 - /j,(t))(x 1 K(x 1 )) *m (a:i) 

m^xijT) = K{xi) * m 2 (xi, r) - 2(1 - [i(t))(xiK(xi)) * mi(ii, r) 

+(1 - ^(r)) 2 ^ 2 ^^!)) * m (a;i) + s 2 K(a;i) * m (xi). 



(Bl) 
(B2) 

(B3) 

(B4a) 
(B4b) 
(B4c) 

(B5a) 
(B5b) 
(B5c) 

(B6a) 
(B6b) 

(B6c) 



Appendix C: Error in the estimation of a 



with an expectation value of a due to Eq. (C2 1 



< V > = a 1 



Let £ be a Gaussian random variable with mean zero and 
standard deviation a. Further let & , . . . £ w be a sample of TV xherefore y can be written 
independent realizations of £. Then the expectation values of 



CiCj an d £ 2 £ 2 ar e g iven b Y 
< fi > = 

< > = cr 2 ^j 

= <r A + 2a% j . 
The variance of £ can be estimated from the sample as 



V 



a 2 + AV 



as 

with 



(C5) 



< AV > = 0. (C6) 



V = 



1 N 



(CI) 
(C2) 
(C3) 



(C4) 



The variance of AV^ can be evaluated using Eq. (C3 i 



<{AVf> = <(^E^-- 2 ) 2 > 

i=l 
N N 



AT 



N 2 



i=l 



i=l j = l 
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1 (iV 2 (T 4 + 2A^cr 4 ) — (T 4 The standard deviation of the error in S is therefore given by 

N 2 

^ (C7) 



Up to now only an estimation of the variance of £ has been 
made. The estimation S of its standard deviation is given by 



{ Avy 



Aa = \ < — > - -= (CIO) 



5 = VV = a,ll + ^ (C8) " ^ V " V™ 



which, for large N, can be approximated as 

1AV, AU 
2^ } = a+ ^a 
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